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Abstract 

We present results of a high resolution numerical study of two dimensional (2d) Rayleigh- Taylor 
turbulence using a recently proposed thermal lattice Boltzmann method (LBT). The goal of our 
study is both methodological and physical. We assess merits and limitations concerning small- and 
large-scale resolution/accuracy of the adopted integration scheme. We discuss quantitatively the 
requirements needed to keep the method stable and precise enough to simulate stratified and un- 
stratified flows driven by thermal active fluctuations at high Rayleigh and high Reynolds numbers. 
We present data with spatial resolution up to 4096 x 10000 grid points and Rayleigh number up 
to Ra ~ 10^^. The statistical quality of the data allows us to investigate velocity and temperature 
fluctuations, scale-by-scale, over roughly four decades. We present a detailed quantitative analysis 
of scaling laws in the viscous, inertial and integral range, supporting the existence of a Bolgiano-like 
inertial scaling, as expected in 2d systems. We also discuss the presence of small/large intermit- 
tent deviation to the scaling of velocity/temperature fluctuations and the Rayleigh dependency of 
gradients flatness. 
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I. INTRODUCTION 



The Rayleigh- Taylor (RT) instability is present whenever we have the superposition of a 



heavy fluid above a lighter one in a constant acceleration field 



EL 



Applications are numerous, 



from inertial-confinement fusion ^ to supernovae explosions 3| and many others Q]. The 
RT instability has been studied for decades, but it still presents several open problems jsl- It 
is important to control the initial and asymptotic evolution of the mixing layer between the 
two miscible fluids; the small-scale turbulent fluctuations, their anisotropic/isotropic ratio; 
their dependency on the initial perturbation spectrum, on the geometry of the containing 
volumes or on the physical dimensions of the embedding space (see [o, l7| for recent high 
resolution numerical studies). Concerning astrophysical and nuclear applications, the two 
fluids evolve with strong compressible and/or stratification effects, a situation which is dif- 
ficult to investigate either theoretically or numerically. The set up studied in this paper is 
two dimensional (2d) and the initial configuration slightly different from what usually found 
in the literature: the spatial temporal evolution of a single component fluid with a cold 
uniform region in the top half and a hot uniform region on the bottom half (see figure [T] 
for details). Such a situation is of interest for convection in the atmosphere, ocean or even 
stars interiors, where, masses of hot /cold fluid may be found in unstable situations js-lOj. 
The choice to focus on a 2d geometry is motivated by different methodological, theoretical 
and phenomenological challenges. First, concerning the method, 2d geometries allow to 
push the numerics to unprecedented resolution - here up to 4096 x 10000 grid points - with 
correspondingly high Rayleigh/ Reynolds numbers; this is an excellent testing ground for the 
lattice Boltzmann Thermal (LBT) scheme [ll,[l2| in fully developed situations, with highly 



intermittent gradient statistics, and a well developed inertial range of scales with power law 
distributions. We initially validate the method against exact relationships originating from 
the hydrodynamical Navier-Stokes- Fourier equations. Then, within the limits settled by the 
validation steps, we show that the scheme -albeit being only second order accurate- allows 
for quantitative studies of hydrodynamical statistical fluctuations over a four decades inter- 
val of scales. This is, to the best of our knowledge, the first time such a huge range of scales 
has ever been explored using LBT codes for turbulent flows. From the phenomenological 
point of view, theoretical work HQ and pioneering numerical simulations at smaller 
resolution tell us that Rayleigh- Taylor dynamics in 2d displays Bolgiano statistics for veloc- 
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ity and temperature fields, at least at scales small enough and far enough from the edges 



of the mixing layer. Bolgiano theory, at variance from Kolmogorov theory [16| , predicts for 
typical inertial-range velocity and temperature fluctuations on a generic inertial scale, i?, 
the following laws 





where L{t) is a measure of the extension of the mixing layer at any given time t during the RT 
evolution, and K{t) is the square root of the total kinetic energy inside the mixing layer (see 
below for a precise definition). These scaling properties tell us that temperature/velocity is 
rougher/smoother than expected for Kolmogorov scaling ~ R^^^. This is due to the active 
role played by buoyancy in the vertical momentum evolution, i.e. temperature becomes 
a fully active scalar at all inertial scales. This is in clear contrast with the Kolmogorov 

13 1 and observed 0, Q in three dimensional (3d) cases. 2d 



like phenomenology expected 



Rayleigh- Taylor systems realize one of those cases where the forcing mechanism -buoyancy- 
overwhelms non-linear energy transfer. This has also theoretical relevance, connected to the 



universality of small-scale statistics in presence of multi-scale forcing mechanisms 181-122 1 



in general, or to renormalization group approaches [23||, in particular. At variance with 
stochastic external multi-scale forcing mechanisms, here the statistics of the buoyancy is 
directly connected to the velocity field itself, opening the way for new phenomena which 
we discuss in details later. Far from being interesting only for theoretical reasons, Bolgiano 
scaling is believed to characterize small scale velocity and temperature fluctuations in 3d 
Rayleigh-Benard convection close to the rigid boundaries, where the viscous and thermal 



boundary layers merge with the bulk region 2^. In fact, thermo-hydrodynamical evolution 



in the proximity of the boundaries is considered to be the key ingredient driving the whole 



cell behavior 



25| 



Here we will be mainly interested to small scale properties, even though large scale evolu- 
tion presents many important open issues, in particular for stratified flows. For example, we 
have recently shown that RT evolution in the set-up of figure [U is stopped by the adiabatic 



gradient in presence of a strongly stratified atmosphere Investigation of small scale 

properties of such situation, as well as the overshooting observed at the edge of the mixing 
layer is in progress and will be reported elsewhere. 



All simulations are performed using an innovative LBT, proposed in 



12| and already val- 



idated concerning large scale properties on the same geometry here investigated 11|. Stable, 



accurate and efficient discrete kinetic methods describing simultaneous hydrodynamical evo- 
lution of momentum and internal energy are notoriously difficult to achieve Q, The 
main difficulties stem from the development of subtle instabilities when the velocity increases 
locally. In recent years, the situation has started to improve, as different attempts have been 
made to describe active thermal modes within a fully discretized Boltzmann approach 



28 



The advantages offered by LB codes are threefold. First, the hydrodynamical manifold 
is described by the whole Navier-Stokes-Fourier equations, with no need to rely on incom- 
pressible and Boussinesq like approximations. Second, the method is particularly efficient 
in dealing with complex bulk or boundary physics, opening the way to incorporate either 
surface tension effects or complex boundary conditions. Last but not least, pressure fluctu- 
ations are fully incorporated in the hydrodynamical evolution, so we do not need to solve 
for Poisson equations; the method becomes fully local in space, allowing for efficient im- 
plementations on massively parallel machines, even if limited interconnection is available. 
Building on this point, our numerical results have been obtained on the QPACE system, 
a massively parallel machine that uses PowerXCell 8i processors connected by a toroidal 



network 
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following the lines of similar older attempts [39 1. 
Results are as follows. In section |TT] we present the notation and the main physical quan- 
tities that we study in this note, including a cursory overview of RT large scale properties. 
In section IIIII we briefly summarize the LBT method, we present the numerical details and 
we discuss the validation steps. In section HV] we present our results on statistical fluctua- 
tions of temperature, velocity, temperature-fluxes and buoyancy terms over the whole range 
of scales accessed by our numerics. We show that velocity statistics is Bolgiano-like with 
very small -if any- intermittent corrections. We discuss the possible origin of these small 
anomalous corrections, in relation with the corresponding small intermittent fluctuations of 
the buoyancy term, a new scenario for 2d turbulence. On the other hand, we show that 
temperature fluctuations are strongly intermittent with high-order moments fully dominated 
by hot /cold fronts. Such strong intermittency has a direct influence also on the tempera- 
ture flux statistics. Our resolution allows us to address quantitatively and scale-by-scale the 



statistical properties of all 
earlier 2d numerical studies 



lydrodynamical fields; this analysis has not been accessible to 



15| and it is still not within reach in the 3d case. Our concluding 
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TABLE I: Parameters for the three types of RT runs. Atwood number At = (T^ — Tu)/{Tfi + T^); 
viscosity u] thermal diffusivity k] gravity g; temperature in the upper half region T^; temperature 
in the lower half region T^; normalization time r = Lx/{g At); adiabatic lenght corresponding to 
the adiabatic gradient = AT/7; dissipative scale calculated at t = r, 77(r); Maximum Rayleigh 
number Ramax'-, number of independent RT evolution Nconf- 



remarks (section |V]) discuss possible further development towards the study of (i) reactive 
Rayleigh- Taylor systems; (ii) strongly stratified systems; (iii) multiphase/multi-component 
Rayleigh- Taylor or convection systems. 



II. RAYLEIGH- TAYLOR SYSTEMS 



The spatio-temporal evolution of a stratified compressible flow, in a external gravity field, 
g > 0, is ruled by the Navier-Stokes- Fourier equations (double indexes are meant summed 
upon) : 

Dtp = -pdiUi 

pDtUi = -diP - pg6i^z + p^djjUi (2) 
pCpDtT - DtP = xdiiT, 
where Dt is the material derivative, /i, x the molecular viscosity and thermal conductivity, 
Cp the specific heat at constant pressure and p, T, P, u are density, temperature, pressure 
and velocity field. Under the assumption that compressibility and stratification are small 
(the situation addressed in this note) and that fluid parameters depend weakly on the local 
thermodynamic fields, one can expand pressure around its hydrostatic value P = Pq + p, 



with dzPo = —gp and p <^ Pq, and perform a small Mach number expansion 

DtUi 
DtT 
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48|: 



diP 
P 



(3) 



kd,{T. 
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FIG. 1: Bottom: Initial configuration for the stratified Rayleigh- Taylor systems. Temperature in 
the upper half is chosen constant T(){z) = Tup while density follow an hydrostatic profile, Pq{z) = 
Pupexp{—g{z — Zc)/Tup), with Zc the central location in the box. In the lower half we have: 
To{z) = Tdown, and po(^) = Pdownexp{-g{z - Zc)/Tdown)- To be at equilibrium, we require to have 
the same pressure at the interface, PupTup = PdownTdown- The temperature jump at the interface 
is smoothed by a tanh profile with a width of the order of 10 grid points. The bold and tiny 
solid lines represent the temperature and density profiles respectively. Top: Snapshot of the RT 
evolution at three times t = (0.5, l,4)r. 

In this approximation, only temperature fluctuations 9 force the system; we have introduced 
the mean temperature T^, kinematic viscosity v = /i/p, thermal diffusivity k = x/i^pP) 
and adiabatic gradient for an ideal gas, 7 = g/cp. The small Mach expansion and small 
stratification decouple the pressure from the internal energy equation, i.e. p in ([3]) is just 
a Lagrange multiplier used to enforce diUi = everywhere. As we will show in the next 
section, the LBT algorithm we are going to use is meant to reproduce the set of equations 
([2]) and ([3]) in the corresponding limit. 

If the adiabatic gradient is negligible, 7 ~ 0, it is well known that starting from an 
unstable initial condition as depicted in figure [1], any small perturbation will lead to a 
turbulent mixture between the hot and cold regions, expanding along the vertical z-direction. 
Concerning large scale quantities, a huge amount of earlier work (e.g. see Ref. [sl) has 
focused on the extimation of the growth rate of the mixing layer extension, L{t), and of the 
total turbulent kinetic energy, K'^it) = -^77^ / dxdzv?, produced by the conversion of the 
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initial potential energy. Using dimensional analysis and self-similar assumptions 
predicts: 

where to is the typical time needed for the system to reach a fully non-linear evo 



40 



4l| one 
(4) 



va. 



ues of the coefficients, a, f3, have been extensively studied both in 2d and 3d [5 



ution. 
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'he 



42l446l| ■ They depend on the definition of L{t), typically taken either as the region where the 



mean temperature profile, averaged over the horizontal direction, T{z) = j dxT{x, z,t), 
is within a given range, for example: T{z) G 0.95[Tup : Tdown], or as an integral property 
over the whole temperature distribution: 

m = f / dxdz e [ ^^^'"'^^"^"^ 1 , (5) 

-^x J L down -'■up 

with e[x] = 2x; < x < 1/2 and e[x] = 2 (1 - x); 1/2 < x < 1. Using the estimate (g]) 
one may predict the whole profile evolution, adopting either simple constant eddy viscosity 
models or more refined Prandtl mixing length theory j49| . In figure [2] we show the growth 
rate of the mixing layer, kinetic energy and the temporal evolution of the temperature profile, 
as an example of typical evolutions of large scale quantities in our numerics. The agreement 
with the expected phenomenology is very satisfactory. Notice a systematic small deviation 
at large times. This deviation is probably due to a transition induced by the evolving aspect 
ratio. When the aspect ratio becomes order one, important horizontal fiuctuations develop 
in the system, preventing an efficient conversion of potential energy in vertical kinetic energy 
(inset of the same figure). 

In this paper, on the other hand, we focus on small scales quantities, i.e. velocity, 
temperature and fiuxes statistics scale-by-scale. In particular we focus on the following 
set of structure functions, based on moments of order p of velocity, temperature or mixed 
increments: 



X, z 



Sif{R,t) = {\6nu,n, 
S^^\R,t) = {\Sne\\SRU,\P) 
sl!\R,t) = {[{6n9)'\6nu.\r/'), 
where we define the increment of a generic hydrodynamical field, A{x,z,t) as SfjA 



(6) 
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L{t)IL, 




FIG. 2: (a) Mean temperature profile at four different times during tlie RT evolution. In the inset 
we show the rescaling according to the instantaneous mixing layer length L{t), {T(z/L{t),t) — 
Tup) / {Tdown — Tup)- The profile rescal es p erfectly and in agreement with the cubic shape predicted 



by a Prandtl mixing length theory [17| (solid line), (b) Evolution of mixing layer, L(t), with 
superposed the best parabolic fit (solid line), using the self-similar prediction (jH). Lower inset: 
ratio between horizontal (n^) and vertical (uj.) kinetic energy, calculated in the whole, half or one 
quarter of the mixing layer: a transition around r ~ 4 is clearly visible. Despite of this slowing 
down, the relative scaling of total kinetic energy with respect to the mixing layer length satisfies 
the scaling This is shown in the upper inset where we have K{t) ~ LV2(t). 



A{x + R, z, t) — A{x, z, t) and the average 

((■)) = ^-— -/ dx dz{-) 
is performed on the whole horizontal direction and on a given vertical range inside the mixing 
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layer. In order to minimize non homogeneous contributions, we typically restrict the vertical 
extension of the averaging region to = \L{t), with L{t) estimated according to the volume 
average ([5]). Moreover, in the correlation functions defined above, we only show the results 
for spatial increments along the fully homogeneous horizontal direction, x. Subscript (B) 
and (F) in the third and fourth row of ([6]) denote the correlation functions driving the time 
evolution of the p-th moment of velocity increments (the buoyancy forcing term) and of the 



temperature flux, respectively. Chertkov in 



13| developed a coherent phenomenology for 



small-scales 2d Rayleigh- Taylor systems, on the reasonable assumptions that (i) the mixing 
layer evolution is adiabatically slow compared to small scales fluctuations; (ii) the amount 
of kinetic energy dissipation at small scales is negligible (absence of direct energy cascade in 
2d turbulence); (iii) temperature is efficiently dissipated at small scales (direct temperature 
cascade). These three ingredients lead to a unique possible dimensional prediction, the 
Bolgiano scaling ([1]). In particular, one expects in the inertial range: 



Sf{R,t) ~ J^P(t)(^)C.(p) r,(t) « R « L{t) 
while in the viscous range: 



S^J,u.{R,t) ~ KP{t)CM\<^(p)(-^\p 



■Lit)' 



s'i\R,t)^K^{t){^)<-(p\^y 
sP{R,t)r^K^/'m^Y^'^K^y 



R < r/(t) 



with 



and 



Ceip) 



p 
5' 



Up) 



Cb{p) = (0(1) + Up)). Cf{p) = (0(2) + 0(1)) 



p 



(7) 



(8) 



(9) 



(10) 



Moreover, according to 2d Bolgiano scaling, the dissipative scale increases with time, as 
r](t) ~ t^/®. The two expressions ([3IHD for inertial and viscous ranges are such that they 
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match at the viscous scale, r]{t). The presence of a non stationary evolution makes the 
problem particularly interesting. The above phenomenology has been already investigated 



numerically in jl5|, where a good agreement with Bolgiano scaling for low order veloc- 
ity structure functions and a departure from Bolgiano dimensional scaling for temperature 
structure functions were measured, for the first time. On one hand, the results presented in 



15| clearly indicates the validity of Chertkov's phenomenology, plus the extra complexity 



of anomalous intermittent corrections to the temperature field. On the other hand, due 



to limited spatial resolution, the authors of 15| could not assess statistical properties in a 



quantitative way scale-by-scale because they had scaling over only about a decade. Our data 
add to the above discussion a detailed investigations of inertial, viscous and integral range 
properties covering all together around 4 decades. We confirm and measure the presence of 
large anomalous corrections to the temperature scaling: 

Ce{p)=p/5 + Ae{p). 

We also show that our data cannot exclude the presence of small deviations from Bolgiano 
scaling also for velocity field, a novel observation, never reported before and somehow sur- 
prising for 2d turbulence: 

= 3p/5 + A„(p). 



III. NUMERICAL METHOD gz VALIDATION STEPS 



A. The thermal lattice Boltzmann algorithm 



In this section we recall the essential features of the computational lattice Boltzmann 
method employed in the numerical simulations. A complete analysis, along with extensive 



validation steps, can be found in 



11 



I2I ]. The thermal- kinetic description of a compressible 



gas/fluid with variable density p, local velocity u, internal energy /C, and subject to a local 
body force density g, is given by the following equations: 



dtp + di{pui) = 
dt{puk) + di{Pik) = pQk 
dtIC + \diqi = pgiUi, 



(11) 
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IG. 3: Scheme of the discrete set of velocities, r is the lattice constant whose value is r ~ 1.1969 



22, 



341 ] ■ To recover the correct degree of isotropy for tensors describing thermal fluxes, one needs 
at least 37 speeds in 2d and 105 speeds in 3d. A smaller set of discrete velocities can be used if 
off-grid vectors are allowed [69 1. 



where Pik and Qi are the momentum and energy fluxes, still unclosed at this level of descrip- 



tion. A recent paper [12|] has shown that it is possible to recover exactly equations (fTTI) . 
starting from a suitable discrete version of the Boltzmann equations with self consistent 
local equilibria. The reference scheme is summarized by the following equation set 



fi{x + ciAt,t + At)-fi{x,t) = -— (fiix,t)-ft'\x,t)) , (12) 

TLB V ^ 



where fi{x,t) represents a probability density function to flnd a particle at space-time lo 



cation {x,t) whose velocity q belongs to a discrete set [33|, IS^I- The Ihs of equation ( 1T2|) 
stands for the streaming step of such probability whereas the rhs represents the relaxation 
towards local Maxwellian distribution function //^^^ with characteristic time tlb- 

The macroscopic flelds (density, momentum and temperature) are defined in terms of the 
lattice Boltzmann populations: 

P = X^/i; pu = ^cifi; DpT = ^\ci - uf fl, (13) 

I I I 

with D the space dimensionality. The novelty of the algorithm here employed stems from 
the form of the equilibrium distribution function. Here, it directly depends on the coarse 
grained variables plus a shift from the local body force term: 



fr = fr ( P, ^ + rLB9, T + ™=__±£^^^ ) . (14) 
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The detailed structure of this equihbrium distribution function can be found in jlll. l33l. |34|. 
Lattice discretization also induces corrections terms in the macroscopic evolution of averaged 
quantities: both momentum and temperature must be renormalized by discretization effects 
in order to recover the correct hydrodynamical description from the discretized lattice Boltz- 
mann variables. The first correction to momentum is given by a pre- and post-collisional 
average 



11 
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Mi: 



u^^> = It H g 



and the first non-trivial correction to the temperature field by 12 1 



rj.{H) ^ rj. 



AD 



Using these "renormalized" hydrodynamical fields it is possible to recover, using a Chapman- 



Enskog expansion [ll|,|l2|, the standard thermo-hydrodynamical equations for a compressible 
fluid with energy conservation. Such procedure, applied to the kinetic equations fllip . sets the 
fluxes Pik and qi equal to their hydrodynamical counterpart describing advection, dissipation 
and diffusion. In two dimensions [D = 2), the resulting equations for the hydrodynamical 



fields are those given in equations ([2]) (for the explicit calculation see ll|). 



B. Details of the numerical simulations 

We use a 2d LBT algorithm, with 37 population fields (the so called D2Q37 model), 
moving in the directions shown in figure |3l We have run on the QPACE Supercomputer 
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a novel massively parallel computer, powered by IBM PowerXCell 8i processors 



(an enhanced version of the Cell processor) that supports our algorithm very efficiently 52 |. 
Three different sets of runs have been performed (parameters are summarized in table H]) 
at varying accuracy: (A) a fully resolved high resolution simulation, up to 4096 x 10000 
collocation points with kinematic viscosity and thermal conductivity large enough to ensure 
optimal resolution of velocity and temperature fields even for large order statistics; (B) a 
less resolved high resolution simulation, up to 4096 x 6000 collocation points, with small 
scale transport parameters a factor 2 smaller than in case (A); (C) a even less resolved case 
with the same resolution of (B) and viscosity a factor 5 smaller than (A). Runs (B) and (C) 
make the Rayleigh and Reynolds numbers as large as possible, even though the statistical 
properties for sub- viscous scales will not be as accurate as for set (A). The remarkable result 
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that we are able to present is that the LBT method is able to reproduce large scale and 
inertial range physics correctly even in those cases (e.g. runs (B) and (C)), where very small 
scales are not resolved correctly. A systematic way to validate the accuracy of the method 
and its convergence towards the hydrodynamical manifold of the kinetic equations is to 
benchmark numerical results against exact relationships coming from the hydrodynamical 
Navier-Stokes equations of motions. For example, large and small scale accuracy can be 
checked via the equations for the kinetic energy and the enstrophy of the systems: 

dt\{v^)v = -ey + g{eu^)v , ^ 

(15) 

dt\kw^)v = -e^ + g{d^Ow)v, 
where the two dissipative terms are = u{{diUjY)v and = iy{w'^)v, and with {{■))v 
we mean the average over the whole volume. These two exact relations probe large and 
small scales, respectively. In figure H] we show the percentage difference between left hand 
side and right hand side normalized with the buoyancy term, for the three set of runs of 
table [H Gradients of each field have been calculated either as a centered difference of the 
hydrodynamical variable, or using the lattice definition 

diA{x) Ri '^wiciA{x + ciAt) 



with wi suitable weights [ll| and q the lattice velocities (see figure[3]); we find that the second 
choice gives better agreement. While the energy balance equation is well verified within a few 
percent for all resolutions, the enstrophy balance for run (B) and (C) is not satisfactory. As a 
result, gradient statistics will be measured only using data from run (A). The next question 
concerns the range of scales at which accuracy becomes acceptable also for runs (B) and 
(C). This can be monitored by plotting a sort of normalized "effective gradient" at different 
scales. In figure |5] we show for temperature and vertical velocity the quantities: Su}{R,t) = 
iLit)/vit))U'^Si'}iR,t)/iKit)R/vit)f and (i?,t) = (L(t)/r/(t))?«(2)^f (i?, t)/(/2/r/(t))2 
at different times during the RT evolution. Clearly, even though run (C) does not resolve 
gradients correctly, i.e. the curves do not reach a well developed plateaux for small scales, 
they superpose well with the well resolved run (A) as soon as i? ~ 5?7(t). This result is 
important and makes us confident that the LBT numerics is quantitatively accurate even 
when small scales are not perfectly smooth, i.e. the method provides for a sort of implicit 
large eddy simulation (ILES) with an effective sub-grid dissipation. 
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FIG. 4: Large and small scales validation of the LBT scheme for the three sets of runs (A)-(C) in 
table I. (a) Difference between the Ihs and the rhs of the energy equation in (jlSp . normalized with 
the buoyancy term, (b) The same of (a) but for the enstrophy equation. 



IV. SMALL SCALES STATISTICS 



As the mixing layer evolves, the effective Rayleigh number, characterizing the thermal 
instability inside the layer, grows. In the presence of stratification the expression for the 
?layleigh number is not unique. It is possible to introduce a ^-dependent Rayleigh number 



(,/f(.))Lnt)(^-7) 
""'^'^'^^ (^/P(.)c.)(./P(.)) ^''^ 
where the notation (■) indicates averages over the horizontal direction. We follow here a 
common procedure defining a Rayleigh number based on the middle plane, i.e. Ra{t) = 
Ra{z = 0,t). Notice that the presence of stratification appears also through the adiabatic 
term 7, i.e. any RT mixing of the kind here studied will be stopped sooner or later once an 
adiabatic atmosphere is reached. For the case when the adiabatic term is not important, 
AT/L(t) ^ 7, the ultimate scaling regime predicted by Kraichnan is expected. In this 
regime there is a relationship between the normalized heat flux and the Rayleigh number 



where the Nusselt number [Nu) is defined as the total heat flux inside the mixing layer 
normalized with its conducting value: Nu = {6uz) / {k/ST j L{t)). Other important output 
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FIG. 5: Normalized effective gradients for runs (A) (o) and (C) (□) at two different times along 
the RT evolution: t = (120000,240000) in LBT units. 
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FIG. 6: Nusselt vs Rayleigh, data from run (A) (o) and (C) (□). Inset: Dimensionless dissipative 
anomalies, ty and e^, at changing time during RT evolution. 
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parameters for the system are kinetic energy dissipation and thermal dissipation, e^^ee- In 
2d we expect that the normahzed = e^/{K^{t)/L{t)) and ie = €e/{iATfK{t)/L{t)) 
vanish and go to a constant for large Rayleigh (Reynolds), respectively. This is tantamount 
to predict the existence of a direct cascade of temperature fluctuations and the absence of a 
kinetic energy dissipation anomaly. The monotonic increase of Rayleigh during the mixing 
layer evolution allows for a check of the previous predictions. In figure [6] we show both 
the Nusselt vs. Rayleigh law, confirming for more than 4 decades the observation of the 
ultimate regime and the behaviour of and eg during the RT evolution. We observe the 
tendency towards a constant non vanishing dissipative anomaly for temperature fluctuations, 
while kinetic energy is becoming smaller and smaller at increasing Rayleigh (Reynolds), as 
expected. 



A. Scale-by-scale statistics 

In figure [7] we show a log-log plot of Sl^\R,t) and Su}{R,t), for different orders and 
different times. We also superpose the inertial range scaling predicted by the dimensional 
Bolgiano prediction. Even though on a log-log scaling the global overall agreement between 
data and dimensional Bolgiano scaling is not bad, important deviations can be seen both 
at the crossover between viscous and inertial range, R ^ ri{t) and around the integral scale, 
R ~ Lit). Let us first investigate the viscous-inertial cross-over. There, typical velocity 
fluctuations have to go from a smooth different iable behaviour 5ru ~ i? to Bolgiano scaling 
jump in the scaling property is therefore not too large, and one must 
expect important sub-leading contributions well inside the inertial range coming from the 
viscous scaling. Such sub-leading term may spoil scaling properties even at high Rayleigh 
values. Differently, for temperature, the jump in the scaling properties from viscous to 
inertial is large (from R^ to R^/^). Sub-leading terms cannot play any role. Anyhow, such 
a big change in the scaling properties cannot happen in a too short range of scales: the 
interval of increments with neither a pure viscous nor a pure inertial scaling should be large 
in this case. 

A "conveniently simple transition function" to encompass both viscous and inertial range 



scaling in a single fitting expression is given by the Batchelor parametrization 



54- 



which 
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FIG. 7: Log-log plot of velocity (a) and temperature (b)scaling for p = 2, 4, 6 at a late time during 
RT evolution (r ~ 5 and data from run (A)). We also plot the corresponding Bolgiano and viscous 
scaling. 



for a generic structure function of order p reads: 

= a — TTT. (17) 

where Cp is a suitable dimensional normalization parameter and Ap is a dimensionless pa- 
rameter taking into account some small possible dependency of the viscous cutoff on the 



order of the correlation function [59H62|. The above expression is the simplest way to glue 
smoothly a differential behaviour in the viscous range, ~ i?^, for R <^ ri{t) with a rough 
scaling, ~ W^^^^ in the inertial range, 'f]{t) <^ R. We need also to match the inertial- integral 
layer, R ~ L(t), where structure functions start to saturate because all hydrodynamical 
fields decorrelate for R ^ L[t). It is easy to generalize the Batchelor parametrization to 

18 



also encompass such a region of scales, reaching a global phenomenological description of 
structure functions valid for all scales: 

= C, — ^ {R" + 5L"(t))-«^)/», (18) 

+ Apffit)) 2 

where in the above expression, the crossover around R ~ L{t) is fixed by the parameter a and 
B (in the following always chosen a = 4, i? = 1). The potentialities of the parametrization 
(|T8|) cannot be appreciated on log-log plots: a detailed scale-by-scale analysis of structure 
functions behaviour is needed. 

A scale-by-scale analysis can be obtained by looking at the so-called Local Scaling Exponents 
(LSE), i.e. the log derivatives of any structure function: 

cilog(i?) ■ ^^^^ 
Whenever we have a pure power law behaviour, the output must be a constant as a function 
of the separation scale, i?, C{j)\R,t) ~ C,{j>,t). The advantage to measure (fT9|) stems from 
the possibility to follow also the cross-over between viscous and inertial range and between 
inertial and integral range, scale-by- scale, hence the name. In figure Owe show for p = 4, 6 
the velocity structure functions against the Batchelor parametrization ( |T8l) for two different 
Rayleigh numbers. In the body of the figure, we plot the LSE for Su]{R,t) from our data 
and superposed with the corresponding expression coming from the parametrization (fTSjl . 
where we have used the Bolgiano value C,u{p) = |p- The agreement is strikingly good; 
considering together all data at different resolutions, we are able to reproduce the viscous, 
inertial and integral scale behaviour over 4 decades of scaling range. The agreement between 
the Bolgiano dimensional prediction and the velocity scaling is very accurate within error 
bars. Notice that the use of LSE with respect to log-log scaling as depicted in the inset of the 
same figure allows to move the discussion from global fit over many orders of magnitude (for 
the latter) to a scale-by- scale fit of 0{1) quantities (for the former). Moving to temperature 
scaling, the scenario changes. In figure M we show the same as figure [8] but for temperature 
and up to p = 8. Here, the agreement with the Batchelor parametrization with the Bolgiano 
dimensional scaling for temperature Ce(p) = is less good, almost acceptable for low order 
moments, but definitely not on top of the numerical data for high order moments. In order to 
achieve a good fit, on the whole range of scale, one needs to introduce anomalous corrections 
to the exponents Ce{p) = p/5 + Ag{p) used in the Batchelor formula. In the same figure we 
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FIG. 8: Local scaling properties for velocity structure functions, Cu^{p\R,t) for run (A) (o) at 
i = 4r and (C) (□) at t = 4r. We show the case with p = 4 (panel a) and p = 6 (panel b). 
Solid lines correspond to the local scaling exponents as predicted from (|18p using the Bolgiano 
dimensional scaling Q, also drawn as an horizontal line of value 12/5 and 18/5 respectively. Error 
bars are calculates out of the scattering between the Nconf different RT evolutions for each run. 
Insets: structure functions, sl}']{R,t), for p = 4 and p = 6 and the two runs (A) and (C) (same 
symbols). The solid line is the parametrization ()18p . 



(eip) 


Bolgiano Ref.[15] here 


p=4 


0.8 0.6 0.6 ± 0.06 


p=6 


1.2 0.7 0.7 ± 0.07 


p=8 


1.6 — 0.8 ± 0.1 



TABLE II: Summary of temperature scaling exponents using the best fit obtained by the 
parametrization (jl8p . using for rj^t) and L{t) the actual values measured on the data 



show indeed how the use of Ae(4) = —0.2, Ag{Q) = —0.5, A5i(8) = —0.8 gives a much better 
agreement between numerical data and the phenomenological parametrization formula ([T8|) . 
This is, in our view, a very clean demonstration of the existence of anomalous scaling for 
temperature fluctuations in 2d RT. The values measured for Ag{p) are in agreement with 
the one presented in 15[. Result on temperature scaling are summarized in table ( IITI) . 

An important feature of RT in 2d, is the active role played by buoyancy at all scales, 
as witnessed by the Bolgiano phenomenology. The interesting point here is that, as the 



20 




- 




Bolgiano 

anomalous 




- 


10-° 
10-1° 


-s'/'\bj) ' ' - 






10-" 


/ RMt) _ 




6/5 




10" 10 10^ 10' 10* 



















10 



10^ 




Rh(t) 



FIG. 9: Same data as in figure [8] but for temperature scaling. We show p = 4, p = 6 and p = 8 
(panels from top to bottom). The solid line corresponds to the parametrization (jlSp using the 
dimensional Bolgiano scaling exponents. The thick solid line is the same parametrization but with 
anomalous scaling exponents. Already for p = 4 and more importantly for p = 6 and p = 8 the LSE 
for the parametrization (jlSp with dimensional Bolgiano scaling (^{p) = p/5 does not fit the numerical 
data. As a guide to the eyes, we also show the Bolgiano inertial range values as horizontal lines in 
each panel. The curves supporting anomalous scaling are obtained with the following correction 
to the exponents: (4) = —0.2, Ag{6) = —0.5 and Ag{8) = —0.7. Insets: structure functions with 
superposed the Batchelor parametrisation with anomalous inertial exponents (solid line). 



buoyancy is driven by temperature fluctuations, the forcing mechanism in the momentum 
equations is given by a non self-similar -intermittent- field. Navier-Stokes equations forced 
with power law forcing, have attracted the attention in the past both for application of 



the renormalization group 23| and for issues concerning small-scales universality, i.e. un- 



derstanding how strong must be the forcing mechanism in order to change the small-scale 
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FIG. 10: Local Scaling Exponent for buoyancy terms, S^^\R,t), with p = 1,3 for run (A) (o) and 
(C) (□). The soUd Une corresponds to the dimensional estimate (jlSp with Cb{p) = Ce(l) + Cuz{p)- 
The two horizontal lines give the expected Bolgiano scaling in the inertial range, 10/5 {p = 3) and 
16/5 {p = 5). 



statistics in turbulent flows |18|-|21|. Typically, for any given system, there exists a critical 
exponent, h^-, characterizing the power law decaying of the forcing spectrum, E{k) ~ k~^, 
such that for b < be the forcing is the leading mechanism of energy exchange at all scales. In 
our case, the very existence of Bolgiano scaling tells us that we fall in the latter class. The 
main interesting differences here, with previous theoretical and numerical studies, is that the 
forcing mechanism is also intermittent, i.e. very different from the typ ical scaling-invariant 



JJ, 



23j. Indeed, the 



Gaussian and delta-correlated in time power-law forcing used in isi 
high intermittency of the temperature scaling shown in the previous section, suggests the 
possibility that some degree of intermittency is also hidden in the velocity field, even though 
a direct measure as the one shown in figure [H] rules out big effects. In figure fTOl we are looking 
directly at the forcing statistics entering in the equation of high order velocity moments, 
what we call the buoyancy structure functions in ([6]), S^^\r, t). As one can see, even there it 
is hard to disentangle any deviations from Bolgiano dimensional scaling. A different scenario 
appears for the temperature flux structure functions, Sp\R,t), as defined in ([6]), shown in 
figure [m Here, a deviation from the dimensional scaling is visible, due to the higher order of 
temperature fields with respect to the velocity fields entering in these correlation functions. 
A possible way to highlight even better intermittent correction is to look at the behaviour 
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FIG. 11: Run (A). Local Scaling Exponent for temperature flux moments, Sp \R, t), with p = 3, 6 
and p = 8 (inset). The solid thin line for p = 3, 6 (main plot) corresponds to the dimensional 
Bolgiano estimate (fTSl) with C,{p) = (C9(2) + Cu2(2))p/3. In the inset we show fit wioth both 
dimensional Bolgiano and anomalous estimates. Notice the better agreement with the anomalous 
case 




of velocity and temperature hyper-flatness: 

Any systematic dependence of flatness on the reference scale R is the signature of a non 
perfect self-similar statistics. Figure [12] shows temperature and velocity flatness at two 
different times during the RT evolution. Temperature is clearly intermittent with a flatness 
which increases at decreasing scale. Velocity is more noisy, nevertheless, our data cannot 
exclude a small scale-dependency of flatness also for the latter, pointing towards small but 
detectable breaking of self-similarity, i.e. corrections to the Bolgiano scaling also for velocity. 
In the inset of the same figure, we show the relative scaling of 4th and 6th order structure 
functions versus the second order one, a procedure known as ESS in literature 
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E4: 



sf\R,t) vs si'\R,ty, Sjf]{R,t) vs 5£)(i?,t); 

Here, a breaking of self similarity is detected as a deviation from the dimensional scaling 
S^^\R,t) = {S^'^\R,t)y^'^. Deviations for the temperature/velocity are strong/small and 
clearly detectable. 
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The above results suggest that in order to highhght some possible non trivial scaling 
properties in the velocity statistics, one needs to look at small scales, where temperature 
intermittency becomes more intense and possibly affects also the momentum equations. In 
figure [T5] we show the behaviour of the flatness of velocity and temperature derivatives during 
the RT evolution, i.e. at increasing Rayleigh: 

p ... m^^y) . p ... mey) .... 

Both small-scales temperature and velocity intermittency are increasing, with the tempera- 
ture case much faster. We fit a power law behaviour: 

F^P) ~ i?a«fi-"-(P); F^^^ ~ Ra^^-'^P^; (22) 

with C,d^u,{p) = 0.12(5) and ^d^e{p) = 0.15(5). While the result for temperature does not 
surprise, the result for velocity does, supporting the existence of a small, but detectable 
intermittent correction to the 2d Bolgiano scaling for the velocity field. 

V. CONCLUSIONS AND FURTHER DEVELOPMENTS 

In this paper we have presented the results of a high resolution numerical study of 2d 
Rayleigh- Taylor turbulence using a new thermal lattice Boltzmann method. The goal of the 
study was both methodological and physical. Concerning the method, we validate and assess 
the stability, accuracy and performances of the numerical discrete kinetic algorithm used, 
showing that even when not perfectly resolved at small scales, inertial and integral scale 
hydrodynamics is well reproduced. This result opens the way to a systematic exploitation of 
LBT algorithms also for fully developed turbulence. Concerning the physics of RT turbulence 
in 2d, we have analyzed data up to Ra ~ 10^^ and shown that the dynamics is dominated by a 
Bolgiano phenomenology, i.e. thermal fluctuations in the buoyancy term are overwhelming 
the kinetic energy flux at all scales. We have also shown that:(i) a suitable Batchelor- 
like parametrization is able to reproduce scale-by-scale the whole statistics at all scales, 
over about 4 decades; (ii) temperature fluctuations show small-scales intermittency, with 
scaling exponents tending to saturate at high orders (see table HIl), a signature of persistence 



of hot/cold fronts even at very small scales 65|; (iii) velocity statistics is much closer to 



Bolgiano dimensional scaling even if small intermittent corrections cannot be ruled out, 
especially concerning gradients evolution. 
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FIG. 12: Velocity (bottom) and temperature (top) flatness at three different times during RT 
evolution t = (2, 3, 5)r in run (A). Insets: ESS plot for structure functions of order p = 6 (+) and 
p = 8 (x) versus structure function with p = 2. The solid line corresponds to the dimensional 
scaling . 



All these results are relevant for 3d thermal systems in presence of boundaries. Indeed, 



Bolgiano physics is believed to describe also thermal and ye, 
boundary in real 3d convective Rayleigh-Benard cells 
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ocity fluctuations close to the 



71| . The existence of anomalous 



intermittent small scale fluctuations also in these cases is relevant to control the physics of 
the viscous and thermal boundary layers. 

The algorithm presented here opens the way for natural generalization to more complex 
situtations. First, it is trivially extendable to 3d cases. Second, it can be further generalized 
including bulk forcing terms in the internal energy equation, to describe reactive system 66 1. 
Third, it is under investigation the possibility to couple the thermal LBT scheme with multi- 
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FIG. 13: Flatness based on velocity □, and temperature A, gradients as a function of the Rayleigh 
number. Two power laws with the best fit for Ra > 10^ are also shown as a guide for the eyes. 



component and/or multi-phase LBT models 67|, including non-trivial wettability properties 



at the boundaries 



681]: a case of interest to describe convection of boiling systems. 
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